rm(list = ls())
require(mfx)
require(zTree)
require(stargazer)
require(stats4)
require(lmtest)
require(bbmle)
require(msm)
require(nlWaldTest)
require(sandwich)


zTT1= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October20ChangingPriors.xls")
#zTT1= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October20ChangingPriors.xls")
subjects1=zTT1$subjects

zTT2= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
#zTT2= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
subjects2=zTT2$subjects
subjects2$Subject=subjects2$Subject+100

zTT3= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October24ChangingPriors.xls")
#zTT3= zTreeTables("C:/Users/nneligh/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October23ChangingPriors.xls")
subjects3=zTT3$subjects
subjects3$Subject=subjects3$Subject+200

subjects=rbind(subjects1,subjects2,subjects3)

priors=c(0.5,0.85)

#Clean data
subjects$responder=(subjects$"response1[50]">0)
subjects=subset(subjects,responder>0)

uidlist=unique(subjects$Subject)

#find N
#**
length(uidlist)
#**


#create a frame with observations
obsframe=data.frame(uid=vector('numeric'),block=vector('numeric'), priorA=vector('numeric'), number=vector('numeric'), state=vector('numeric'), response=vector('numeric'))
for(i in 1:length(uidlist)){
playframe=subjects[i,]
uidplace=uidlist[i]

#block 1
for (j in 1:50){
blockplace=1
priorplace=priors[playframe[1,682]]
numberplace=j
stateplace=playframe[1,(735+j)]
responseplace=playframe[1,(71+j)]
addframe=data.frame(uid=uidplace,block=blockplace, priorA=priorplace, number=numberplace, state=stateplace, response=responseplace)
obsframe=rbind(obsframe,addframe)
}

#block 2
for (j in 1:50){i
blockplace=2
priorplace=priors[playframe[1,683]]
numberplace=j
stateplace=playframe[1,(785+j)]
responseplace=playframe[1,(221+j)]
addframe=data.frame(uid=uidplace,block=blockplace, priorA=priorplace, number=numberplace, state=stateplace, response=responseplace)
obsframe=rbind(obsframe,addframe)

}

}

expframe=obsframe

#definitions
expframe$priorB=1-expframe$priorA
expframe$Aresponse=(expframe$response==1)
expframe$Bresponse=(expframe$response==2)
expframe$stateA=(expframe$state==1)
expframe$stateB=(expframe$state==2)
expframe$correct=(expframe$response==expframe$state)

write.table(expframe, "C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Revision/Reports/AmbhuelChangingPriors.csv", sep = ",", col.names = T, append = F)


#fancy regression for easy clustering by Markus Conrad https://www.r-bloggers.com/2021/05/clustered-standard-errors-with-r/
ols <- function(form, data, robust=FALSE, cluster=NULL,digits=3){
  r1 <- lm(form, data)
  if(length(cluster)!=0){
    data <- na.omit(data[,c(colnames(r1$model),cluster)])
    r1 <- lm(form, data)
  }
  X <- model.matrix(r1)
  n <- dim(X)[1]
  k <- dim(X)[2]
  if(robust==FALSE & length(cluster)==0){
    se <- sqrt(diag(solve(crossprod(X)) * as.numeric(crossprod(resid(r1))/(n-k))))
    res <- cbind(coef(r1),se)
  }
  if(robust==TRUE){
    u <- matrix(resid(r1))
    meat1 <- t(X) %*% diag(diag(crossprod(t(u)))) %*% X
    dfc <- n/(n-k)    
    se <- sqrt(dfc*diag(solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X))))
    vcov <- solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X))
    res <- cbind(coef(r1),se)
    }
  if(length(cluster)!=0){
    clus <- cbind(X,data[,cluster],resid(r1))
    colnames(clus)[(dim(clus)[2]-1):dim(clus)[2]] <- c(cluster,"resid")
    m <- dim(table(clus[,cluster]))
    dfc <- (m/(m-1))*((n-1)/(n-k))
    uclust  <- apply(resid(r1)*X,2, function(x) tapply(x, clus[,cluster], sum))
    se <- sqrt(diag(solve(crossprod(X)) %*% (t(uclust) %*% uclust) %*% solve(crossprod(X)))*dfc)   
    vcov <- solve(crossprod(X)) %*% (t(uclust) %*% uclust) %*% solve(crossprod(X))
    res <- cbind(coef(r1),se)
    }
  res <- cbind(res,res[,1]/res[,2],(1-pnorm(abs(res[,1]/res[,2])))*2)
  res1 <- matrix(as.numeric(sprintf(paste("%.",paste(digits,"f",sep=""),sep=""),res)),nrow=dim(res)[1])
  rownames(res1) <- rownames(res)
  colnames(res1) <- c("Estimate","Std. Error","t value","Pr(>|t|)")
  out=list(res1,vcov)
  return(out)
}


#Find thresholds and look for dropping

priorchangeframe=data.frame(uid=vector('numeric'),pA0.5=vector('numeric'),pA0.5lower=vector('numeric'),pA0.5upper=vector('numeric'),pA0.85=vector('numeric'),ChangePost0.5_0.85=vector('numeric'),Chisq0.5_0.85=vector('numeric'),Prob0.5_0.85=vector('numeric'))

for(i in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[i])
model=lm(Aresponse~factor(stateA+priorA)+0,data=playframe)
na1=model$coef[1]
na2=model$coef[2]
a1=model$coef[3]
a2=model$coef[4]


if(sum(playframe$Aresponse*(playframe$priorA==0.85))==0){
postAhigh0.85=0
}else{
postAhigh0.85=(0.85*a2)/(0.85*a2+0.15*na2)
}

if(sum(playframe$Aresponse*(playframe$priorA==0.5))==0){

postAhigh0.5=0
postAhigh0.5upper=0
postAhigh0.5lower=0

diff0.5_0.85=postAhigh0.5-postAhigh0.85
chisq0.5_0.85=1e7*(postAhigh0.5-postAhigh0.85>0)
diff0.5_0.85p=7.888609e-31*(postAhigh0.5-postAhigh0.85>0)+(1-7.888609e-31)*(postAhigh0.5-postAhigh0.85==0)


}else{
postAhigh0.5=(0.5*a1)/(0.5*a1+0.5*na1)

vcovmod=vcovHC(model, type="HC3")

threshsterr=as.numeric(deltamethod(~(0.5*x1)/(0.5*x1+0.5*x2),coef(model)[c(3,1)],vcovmod[c(3,1),c(3,1)]))
postAhigh0.5upper=postAhigh0.5+1.96*threshsterr
postAhigh0.5lower=postAhigh0.5-1.96*threshsterr

if(sum(playframe$Aresponse)==length(playframe$Aresponse)|sum(playframe$Aresponse)==0){

diff0.5_0.85=0
chisq0.5_0.85=0
diff0.5_0.85p=1

}else{

diff0.5_0.85=postAhigh0.5-postAhigh0.85

if(diff0.5_0.85==0){
chisq0.5_0.85=0
diff0.5_0.85p=1
}else{
test0.5_0.85=nlWaldtest(model,"(0.5*b[3])/(0.5*b[3]+0.5*b[1])-(0.85*b[4])/(0.85*b[4]+0.15*b[2])",Vcov=vcovmod)
chisq0.5_0.85=as.numeric(test0.5_0.85[1])
diff0.5_0.85p=as.numeric(test0.5_0.85[4])
}
}

}
addframe=data.frame(uid=uidlist[i],pA0.5=postAhigh0.5,pA0.5lower=postAhigh0.5lower,pA0.5upper=postAhigh0.5upper,pA0.85=postAhigh0.85,ChangePost0.5_0.85=diff0.5_0.85,Chisq0.5_0.85=chisq0.5_0.85,Prob0.5_0.85=diff0.5_0.85p)


priorchangeframe=rbind(priorchangeframe,addframe)
}

percent=vector('numeric')
percent[1]=100*sum((priorchangeframe$pA0.5<0.85))/length(priorchangeframe$uid)
percent[2]=100*sum((priorchangeframe$pA0.5>=0.85)*(priorchangeframe$pA0.5<=1))/length(priorchangeframe$uid)
thresholds=c("[0.5,0.85)","[0.85,1]")

thresholdsframe=data.frame(Theshold=thresholds,Percent=percent)

#Make some graphs of posteriors with subjects that have different thresholds

#Next Graph
above0.85=subset(priorchangeframe, pA0.5lower>0.85)
uidlistabove0.85=unique(above0.85$uid)
above0.85frame0.5=subset(expframe,priorA==0.5&uid %in%uidlistabove0.85)
above0.85frame0.85=subset(expframe,priorA==0.85&uid %in%uidlistabove0.85)

#for p values

comboframe0.85=rbind(above0.85frame0.5,above0.85frame0.85)
comboframe0.85$highpri=(comboframe0.85$priorA==0.85)
diffmodelA=ols(stateB~Aresponse+highpri+highpri*Aresponse,data=comboframe0.85,robust=TRUE, cluster="uid")
diffmodelB=ols(stateA~Bresponse+highpri+highpri*Bresponse,data=comboframe0.85,robust=TRUE, cluster="uid")

#number of significant violations
priorchangeframeabove0.85=subset(priorchangeframe,uid %in% uidlistabove0.85)

#****************************************
sum(priorchangeframeabove0.85$Prob0.5_0.85<0.05)/length(priorchangeframeabove0.85$Prob0.5_0.85)
#****************************************

#for graphs

model0.85B0.5=ols(stateB~Aresponse, data=above0.85frame0.5, robust=TRUE, cluster="uid")[[1]]
model0.85A0.5=ols(stateA~Bresponse, data=above0.85frame0.5, robust=TRUE, cluster="uid")[[1]]
model0.85B0.85=ols(stateB~Aresponse, data=above0.85frame0.85, robust=TRUE, cluster="uid")[[1]]
model0.85A0.85=ols(stateA~Bresponse, data=above0.85frame0.85, robust=TRUE, cluster="uid")[[1]]


correctfrac0.5A=sum(above0.85frame0.5$correct*above0.85frame0.5$Aresponse)/sum(above0.85frame0.5$Aresponse)
correctfrac0.5B=sum(above0.85frame0.5$correct*above0.85frame0.5$Bresponse)/sum(above0.85frame0.5$Bresponse)
correctfrac0.85A=sum(above0.85frame0.85$correct*above0.85frame0.85$Aresponse)/sum(above0.85frame0.85$Aresponse)
correctfrac0.85B=sum(above0.85frame0.85$correct*above0.85frame0.85$Bresponse)/sum(above0.85frame0.85$Bresponse)

correctfrac0.5=c(correctfrac0.5A,correctfrac0.5B)
correctfrac0.85=c(correctfrac0.85A,correctfrac0.85B)
correctfrac0.5high=correctfrac0.5+c(model0.85A0.5[1,2],model0.85B0.5[1,2])
correctfrac0.85high=correctfrac0.85+c(model0.85A0.85[1,2],model0.85B0.85[1,2])
correctfrac0.5low=correctfrac0.5-c(model0.85A0.5[1,2],model0.85B0.5[1,2])
correctfrac0.85low=correctfrac0.85-c(model0.85A0.85[1,2],model0.85B0.85[1,2])

#*****************************
barCenters=barplot(matrix(c(rbind(correctfrac0.5,correctfrac0.85)),nr=2), beside=T, 
        col=c("blue","orange"), 
        names.arg=c("P(1|A)","P(2|B)"),ylim = c(0, 1),xpd = FALSE)
 par(xpd=TRUE)
legend(2.75,-0.05,ncol=2, title="Prior", c("0.5","0.85"), pch=15, 
       col=c("blue","orange"), 
       bty="n")

box(bty="l")

segments(barCenters,c(rbind(correctfrac0.5low,correctfrac0.85low)), barCenters,
        c(rbind(correctfrac0.5high,correctfrac0.85high)), lwd = 1.5)

arrows(barCenters, c(rbind(correctfrac0.5low,correctfrac0.85low)), barCenters,
        c(rbind(correctfrac0.5high,correctfrac0.85high)), lwd = 1.5, angle = 90,
       code = 3, length = 0.05)

#*******************************


#report p-values for differences
above0.85frame=subset(expframe,uid %in% uidlistabove0.85 & (priorA==0.5|priorA==0.85))
above0.85frame$priorA0.85=(above0.85frame$priorA==0.85)
model0.85B=ols(stateB~Aresponse+priorA0.85+priorA0.85*Aresponse, data=above0.85frame, robust=TRUE, cluster="uid")[[1]]
model0.85A=ols(stateA~Bresponse+priorA0.85+priorA0.85*Bresponse, data=above0.85frame, robust=TRUE, cluster="uid")[[1]]

#**********************
model0.85B[3,1]
model0.85B[3,4]
model0.85A[3,1]
model0.85A[3,4]
#**********************

